<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>zsteqr.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="ZSTEQR.1"></a><a href="zsteqr.f.html#ZSTEQR.1">ZSTEQR</a>( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          COMPZ
      INTEGER            INFO, LDZ, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      DOUBLE PRECISION   D( * ), E( * ), WORK( * )
      COMPLEX*16         Z( LDZ, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="ZSTEQR.19"></a><a href="zsteqr.f.html#ZSTEQR.1">ZSTEQR</a> computes all eigenvalues and, optionally, eigenvectors of a
</span><span class="comment">*</span><span class="comment">  symmetric tridiagonal matrix using the implicit QL or QR method.
</span><span class="comment">*</span><span class="comment">  The eigenvectors of a full or band complex Hermitian matrix can also
</span><span class="comment">*</span><span class="comment">  be found if <a name="ZHETRD.22"></a><a href="zhetrd.f.html#ZHETRD.1">ZHETRD</a> or <a name="ZHPTRD.22"></a><a href="zhptrd.f.html#ZHPTRD.1">ZHPTRD</a> or <a name="ZHBTRD.22"></a><a href="zhbtrd.f.html#ZHBTRD.1">ZHBTRD</a> has been used to reduce this
</span><span class="comment">*</span><span class="comment">  matrix to tridiagonal form.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  COMPZ   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'N':  Compute eigenvalues only.
</span><span class="comment">*</span><span class="comment">          = 'V':  Compute eigenvalues and eigenvectors of the original
</span><span class="comment">*</span><span class="comment">                  Hermitian matrix.  On entry, Z must contain the
</span><span class="comment">*</span><span class="comment">                  unitary matrix used to reduce the original matrix
</span><span class="comment">*</span><span class="comment">                  to tridiagonal form.
</span><span class="comment">*</span><span class="comment">          = 'I':  Compute eigenvalues and eigenvectors of the
</span><span class="comment">*</span><span class="comment">                  tridiagonal matrix.  Z is initialized to the identity
</span><span class="comment">*</span><span class="comment">                  matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the diagonal elements of the tridiagonal matrix.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, the eigenvalues in ascending order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, the (n-1) subdiagonal elements of the tridiagonal
</span><span class="comment">*</span><span class="comment">          matrix.
</span><span class="comment">*</span><span class="comment">          On exit, E has been destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N)
</span><span class="comment">*</span><span class="comment">          On entry, if  COMPZ = 'V', then Z contains the unitary
</span><span class="comment">*</span><span class="comment">          matrix used in the reduction to tridiagonal form.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
</span><span class="comment">*</span><span class="comment">          orthonormal eigenvectors of the original Hermitian matrix,
</span><span class="comment">*</span><span class="comment">          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
</span><span class="comment">*</span><span class="comment">          of the symmetric tridiagonal matrix.
</span><span class="comment">*</span><span class="comment">          If COMPZ = 'N', then Z is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDZ     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array Z.  LDZ &gt;= 1, and if
</span><span class="comment">*</span><span class="comment">          eigenvectors are desired, then  LDZ &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
</span><span class="comment">*</span><span class="comment">          If COMPZ = 'N', then WORK is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  the algorithm has failed to find all the eigenvalues in
</span><span class="comment">*</span><span class="comment">                a total of 30*N iterations; if INFO = i, then i
</span><span class="comment">*</span><span class="comment">                elements of E have not converged to zero; on exit, D
</span><span class="comment">*</span><span class="comment">                and E contain the elements of a symmetric tridiagonal
</span><span class="comment">*</span><span class="comment">                matrix which is unitarily similar to the original
</span><span class="comment">*</span><span class="comment">                matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE, TWO, THREE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                   THREE = 3.0D0 )
      COMPLEX*16         CZERO, CONE
      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
     $                   CONE = ( 1.0D0, 0.0D0 ) )
      INTEGER            MAXIT
      PARAMETER          ( MAXIT = 30 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
     $                   LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
     $                   NM1, NMAXIT
      DOUBLE PRECISION   ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
     $                   S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.96"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.97"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.97"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>, <a name="DLAPY2.97"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
      EXTERNAL           <a name="LSAME.98"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.98"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="DLANST.98"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>, <a name="DLAPY2.98"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLAE2.101"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>, <a name="DLAEV2.101"></a><a href="dlaev2.f.html#DLAEV2.1">DLAEV2</a>, <a name="DLARTG.101"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>, <a name="DLASCL.101"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>, <a name="DLASRT.101"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>, <a name="XERBLA.101"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>,
     $                   <a name="ZLASET.102"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>, <a name="ZLASR.102"></a><a href="zlasr.f.html#ZLASR.1">ZLASR</a>, ZSWAP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, SIGN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.113"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'N'</span> ) ) THEN
         ICOMPZ = 0
      ELSE IF( <a name="LSAME.115"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'V'</span> ) ) THEN
         ICOMPZ = 1
      ELSE IF( <a name="LSAME.117"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPZ, <span class="string">'I'</span> ) ) THEN
         ICOMPZ = 2
      ELSE
         ICOMPZ = -1
      END IF
      IF( ICOMPZ.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
     $         N ) ) ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.131"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZSTEQR.131"></a><a href="zsteqr.f.html#ZSTEQR.1">ZSTEQR</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         IF( ICOMPZ.EQ.2 )
     $      Z( 1, 1 ) = CONE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine the unit roundoff and over/underflow thresholds.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.148"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'E'</span> )
      EPS2 = EPS**2
      SAFMIN = <a name="DLAMCH.150"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SAFMAX = ONE / SAFMIN
      SSFMAX = SQRT( SAFMAX ) / THREE
      SSFMIN = SQRT( SAFMIN ) / EPS2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the eigenvalues and eigenvectors of the tridiagonal
</span><span class="comment">*</span><span class="comment">     matrix.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ICOMPZ.EQ.2 )
     $   CALL <a name="ZLASET.159"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, N, N, CZERO, CONE, Z, LDZ )
<span class="comment">*</span><span class="comment">
</span>      NMAXIT = N*MAXIT
      JTOT = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine where the matrix splits and choose QL or QR iteration
</span><span class="comment">*</span><span class="comment">     for each block, according to whether top or bottom diagonal
</span><span class="comment">*</span><span class="comment">     element is smaller.
</span><span class="comment">*</span><span class="comment">
</span>      L1 = 1
      NM1 = N - 1
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      IF( L1.GT.N )
     $   GO TO 160
      IF( L1.GT.1 )
     $   E( L1-1 ) = ZERO
      IF( L1.LE.NM1 ) THEN
         DO 20 M = L1, NM1
            TST = ABS( E( M ) )
            IF( TST.EQ.ZERO )
     $         GO TO 30
            IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
     $          1 ) ) ) )*EPS ) THEN
               E( M ) = ZERO
               GO TO 30
            END IF
   20    CONTINUE
      END IF
      M = N
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
      L = L1
      LSV = L
      LEND = M
      LENDSV = LEND
      L1 = M + 1
      IF( LEND.EQ.L )
     $   GO TO 10
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale submatrix in rows and columns L to LEND
</span><span class="comment">*</span><span class="comment">
</span>      ANORM = <a name="DLANST.201"></a><a href="dlanst.f.html#DLANST.1">DLANST</a>( <span class="string">'I'</span>, LEND-L+1, D( L ), E( L ) )
      ISCALE = 0
      IF( ANORM.EQ.ZERO )
     $   GO TO 10
      IF( ANORM.GT.SSFMAX ) THEN
         ISCALE = 1
         CALL <a name="DLASCL.207"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="DLASCL.209"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
     $                INFO )
      ELSE IF( ANORM.LT.SSFMIN ) THEN
         ISCALE = 2
         CALL <a name="DLASCL.213"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
     $                INFO )
         CALL <a name="DLASCL.215"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
     $                INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Choose between QL and QR iteration
</span><span class="comment">*</span><span class="comment">
</span>      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
         LEND = LSV
         L = LENDSV
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( LEND.GT.L ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QL Iteration
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small subdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>   40    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDM1 = LEND - 1
            DO 50 M = L, LENDM1
               TST = ABS( E( M ) )**2
               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M+1 ) )+
     $             SAFMIN )GO TO 60
   50       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         M = LEND
<span class="comment">*</span><span class="comment">
</span>   60    CONTINUE
         IF( M.LT.LEND )
     $      E( M ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 80
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If remaining matrix is 2-by-2, use <a name="DLAE2.251"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a> or <a name="SLAEV2.251"></a><a href="slaev2.f.html#SLAEV2.1">SLAEV2</a>
</span><span class="comment">*</span><span class="comment">        to compute its eigensystem.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L+1 ) THEN
            IF( ICOMPZ.GT.0 ) THEN
               CALL <a name="DLAEV2.256"></a><a href="dlaev2.f.html#DLAEV2.1">DLAEV2</a>( D( L ), E( L ), D( L+1 ), RT1, RT2, C, S )
               WORK( L ) = C
               WORK( N-1+L ) = S
               CALL <a name="ZLASR.259"></a><a href="zlasr.f.html#ZLASR.1">ZLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, 2, WORK( L ),
     $                     WORK( N-1+L ), Z( 1, L ), LDZ )
            ELSE
               CALL <a name="DLAE2.262"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>( D( L ), E( L ), D( L+1 ), RT1, RT2 )
            END IF
            D( L ) = RT1
            D( L+1 ) = RT2
            E( L ) = ZERO
            L = L + 2
            IF( L.LE.LEND )
     $         GO TO 40
            GO TO 140
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 140
         JTOT = JTOT + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Form shift.
</span><span class="comment">*</span><span class="comment">
</span>         G = ( D( L+1 )-P ) / ( TWO*E( L ) )
         R = <a name="DLAPY2.280"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( G, ONE )
         G = D( M ) - P + ( E( L ) / ( G+SIGN( R, G ) ) )
<span class="comment">*</span><span class="comment">
</span>         S = ONE
         C = ONE
         P = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Inner loop
</span><span class="comment">*</span><span class="comment">
</span>         MM1 = M - 1
         DO 70 I = MM1, L, -1
            F = S*E( I )
            B = C*E( I )
            CALL <a name="DLARTG.293"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( G, F, C, S, R )
            IF( I.NE.M-1 )
     $         E( I+1 ) = R
            G = D( I+1 ) - P
            R = ( D( I )-G )*S + TWO*C*B
            P = S*R
            D( I+1 ) = G + P
            G = C*R - B
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If eigenvectors are desired, then save rotations.
</span><span class="comment">*</span><span class="comment">
</span>            IF( ICOMPZ.GT.0 ) THEN
               WORK( I ) = C
               WORK( N-1+I ) = -S
            END IF
<span class="comment">*</span><span class="comment">
</span>   70    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If eigenvectors are desired, then apply saved rotations.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ICOMPZ.GT.0 ) THEN
            MM = M - L + 1
            CALL <a name="ZLASR.315"></a><a href="zlasr.f.html#ZLASR.1">ZLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, N, MM, WORK( L ), WORK( N-1+L ),
     $                  Z( 1, L ), LDZ )
         END IF
<span class="comment">*</span><span class="comment">
</span>         D( L ) = D( L ) - P
         E( L ) = G
         GO TO 40
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvalue found.
</span><span class="comment">*</span><span class="comment">
</span>   80    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L + 1
         IF( L.LE.LEND )
     $      GO TO 40
         GO TO 140
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        QR Iteration
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Look for small superdiagonal element.
</span><span class="comment">*</span><span class="comment">
</span>   90    CONTINUE
         IF( L.NE.LEND ) THEN
            LENDP1 = LEND + 1
            DO 100 M = L, LENDP1, -1
               TST = ABS( E( M-1 ) )**2
               IF( TST.LE.( EPS2*ABS( D( M ) ) )*ABS( D( M-1 ) )+
     $             SAFMIN )GO TO 110
  100       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         M = LEND
<span class="comment">*</span><span class="comment">
</span>  110    CONTINUE
         IF( M.GT.LEND )
     $      E( M-1 ) = ZERO
         P = D( L )
         IF( M.EQ.L )
     $      GO TO 130
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If remaining matrix is 2-by-2, use <a name="DLAE2.358"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a> or <a name="SLAEV2.358"></a><a href="slaev2.f.html#SLAEV2.1">SLAEV2</a>
</span><span class="comment">*</span><span class="comment">        to compute its eigensystem.
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.EQ.L-1 ) THEN
            IF( ICOMPZ.GT.0 ) THEN
               CALL <a name="DLAEV2.363"></a><a href="dlaev2.f.html#DLAEV2.1">DLAEV2</a>( D( L-1 ), E( L-1 ), D( L ), RT1, RT2, C, S )
               WORK( M ) = C
               WORK( N-1+M ) = S
               CALL <a name="ZLASR.366"></a><a href="zlasr.f.html#ZLASR.1">ZLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, N, 2, WORK( M ),
     $                     WORK( N-1+M ), Z( 1, L-1 ), LDZ )
            ELSE
               CALL <a name="DLAE2.369"></a><a href="dlae2.f.html#DLAE2.1">DLAE2</a>( D( L-1 ), E( L-1 ), D( L ), RT1, RT2 )
            END IF
            D( L-1 ) = RT1
            D( L ) = RT2
            E( L-1 ) = ZERO
            L = L - 2
            IF( L.GE.LEND )
     $         GO TO 90
            GO TO 140
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( JTOT.EQ.NMAXIT )
     $      GO TO 140
         JTOT = JTOT + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Form shift.
</span><span class="comment">*</span><span class="comment">
</span>         G = ( D( L-1 )-P ) / ( TWO*E( L-1 ) )
         R = <a name="DLAPY2.387"></a><a href="dlapy2.f.html#DLAPY2.1">DLAPY2</a>( G, ONE )
         G = D( M ) - P + ( E( L-1 ) / ( G+SIGN( R, G ) ) )
<span class="comment">*</span><span class="comment">
</span>         S = ONE
         C = ONE
         P = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Inner loop
</span><span class="comment">*</span><span class="comment">
</span>         LM1 = L - 1
         DO 120 I = M, LM1
            F = S*E( I )
            B = C*E( I )
            CALL <a name="DLARTG.400"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( G, F, C, S, R )
            IF( I.NE.M )
     $         E( I-1 ) = R
            G = D( I ) - P
            R = ( D( I+1 )-G )*S + TWO*C*B
            P = S*R
            D( I ) = G + P
            G = C*R - B
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If eigenvectors are desired, then save rotations.
</span><span class="comment">*</span><span class="comment">
</span>            IF( ICOMPZ.GT.0 ) THEN
               WORK( I ) = C
               WORK( N-1+I ) = S
            END IF
<span class="comment">*</span><span class="comment">
</span>  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If eigenvectors are desired, then apply saved rotations.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ICOMPZ.GT.0 ) THEN
            MM = L - M + 1
            CALL <a name="ZLASR.422"></a><a href="zlasr.f.html#ZLASR.1">ZLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, N, MM, WORK( M ), WORK( N-1+M ),
     $                  Z( 1, M ), LDZ )
         END IF
<span class="comment">*</span><span class="comment">
</span>         D( L ) = D( L ) - P
         E( LM1 ) = G
         GO TO 90
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Eigenvalue found.
</span><span class="comment">*</span><span class="comment">
</span>  130    CONTINUE
         D( L ) = P
<span class="comment">*</span><span class="comment">
</span>         L = L - 1
         IF( L.GE.LEND )
     $      GO TO 90
         GO TO 140
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling if necessary
</span><span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
      IF( ISCALE.EQ.1 ) THEN
         CALL <a name="DLASCL.446"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
         CALL <a name="DLASCL.448"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMAX, ANORM, LENDSV-LSV, 1, E( LSV ),
     $                N, INFO )
      ELSE IF( ISCALE.EQ.2 ) THEN
         CALL <a name="DLASCL.451"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
     $                D( LSV ), N, INFO )
         CALL <a name="DLASCL.453"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SSFMIN, ANORM, LENDSV-LSV, 1, E( LSV ),
     $                N, INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for no convergence to an eigenvalue after a total
</span><span class="comment">*</span><span class="comment">     of N*MAXIT iterations.
</span><span class="comment">*</span><span class="comment">
</span>      IF( JTOT.EQ.NMAXIT ) THEN
         DO 150 I = 1, N - 1
            IF( E( I ).NE.ZERO )
     $         INFO = INFO + 1
  150    CONTINUE
         RETURN
      END IF
      GO TO 10
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Order eigenvalues and eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      IF( ICOMPZ.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use Quick Sort
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLASRT.476"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>( <span class="string">'I'</span>, N, D, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use Selection Sort to minimize swaps of eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>         DO 180 II = 2, N
            I = II - 1
            K = I
            P = D( I )
            DO 170 J = II, N
               IF( D( J ).LT.P ) THEN
                  K = J
                  P = D( J )
               END IF
  170       CONTINUE
            IF( K.NE.I ) THEN
               D( K ) = D( I )
               D( I ) = P
               CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
            END IF
  180    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZSTEQR.501"></a><a href="zsteqr.f.html#ZSTEQR.1">ZSTEQR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
